home *** CD-ROM | disk | FTP | other *** search
/ Aminet 30 / Aminet 30 (1999)(Schatztruhe)[!][Apr 1999].iso / Aminet / gfx / misc / gnuplot-3.7src.lha / gnuplot-3.7src / gnuplot-3.7.lha / gnuplot-3.7 / contour.c < prev    next >
Encoding:
C/C++ Source or Header  |  1998-12-04  |  45.9 KB  |  1,446 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: contour.c,v 1.31 1998/04/14 00:15:15 drd Exp $";
  3. #endif
  4.  
  5. /* GNUPLOT - contour.c */
  6.  
  7. /*[
  8.  * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
  9.  *
  10.  * Permission to use, copy, and distribute this software and its
  11.  * documentation for any purpose with or without fee is hereby granted,
  12.  * provided that the above copyright notice appear in all copies and
  13.  * that both that copyright notice and this permission notice appear
  14.  * in supporting documentation.
  15.  *
  16.  * Permission to modify the software is granted, but not the right to
  17.  * distribute the complete modified source code.  Modifications are to
  18.  * be distributed as patches to the released version.  Permission to
  19.  * distribute binaries produced by compiling modified sources is granted,
  20.  * provided you
  21.  *   1. distribute the corresponding source modifications from the
  22.  *    released version in the form of a patch file along with the binaries,
  23.  *   2. add special version identification to distinguish your version
  24.  *    in addition to the base release version number,
  25.  *   3. provide your name and address as the primary contact for the
  26.  *    support of your modified version, and
  27.  *   4. retain our contact information in regard to use of the base
  28.  *    software.
  29.  * Permission to distribute the released version of the source code along
  30.  * with corresponding source modifications in the form of a patch file is
  31.  * granted with same provisions 2 through 4 for binary distributions.
  32.  *
  33.  * This software is provided "as is" without express or implied warranty
  34.  * to the extent permitted by applicable law.
  35. ]*/
  36.  
  37.  
  38. /* 
  39.  * AUTHORS
  40.  * 
  41.  *   Original Software:
  42.  *       Gershon Elber
  43.  * 
  44.  *   Improvements to the numerical algorithms:
  45.  *        Hans-Martin Keller, 1995,1997 (hkeller@gwdg.de)
  46.  *
  47.  */
  48.  
  49. #include "plot.h"
  50. #include "setshow.h"
  51.  
  52. #define DEFAULT_NUM_APPROX_PTS 5
  53. #define DEFAULT_BSPLINE_ORDER  3
  54. #define MAX_NUM_APPROX_PTS     100
  55. #define MAX_BSPLINE_ORDER      10    /* ?? Not used ?? */
  56.  
  57. /* for some reason these symbols are also defined in plot.h under different */
  58. /* names */
  59. #define INTERP_NOTHING  CONTOUR_KIND_LINEAR    /* Kind of interpolations on contours. */
  60. #define INTERP_CUBIC    CONTOUR_KIND_CUBIC_SPL    /* Cubic spline interp. */
  61. #define APPROX_BSPLINE  CONTOUR_KIND_BSPLINE    /* Bspline interpolation. */
  62.  
  63. #define ACTIVE     1        /* Status of edges at certain Z level. */
  64. #define INACTIVE   2
  65. #define INNER_MESH 1        /* position of edge in mesh */
  66. #define BOUNDARY   2
  67. #define DIAGONAL   3
  68.  
  69. #define OPEN_CONTOUR     1    /* Contour kinds. */
  70. #define CLOSED_CONTOUR   2
  71.  
  72. #define EPSILON  1e-5        /* Used to decide if two float are equal. */
  73.  
  74. #ifndef TRUE
  75. #define TRUE     -1
  76. #define FALSE    0
  77. #endif
  78.  
  79. #define MAX_POINTS_PER_CNTR     100
  80.  
  81. #define ABS(x)  ((x) > 0 ? (x) : (-(x)))
  82. #define SQR(x)  ((x) * (x))
  83.  
  84. /*
  85.  * struct vrtx_struct {
  86.  *      double X, Y, Z;
  87.  *      struct vrtx_struct *next;
  88.  * };
  89.  * 
  90.  * replaced by 'struct coordinate GPHUGE ', see plot.h (HMK 1997) 
  91.  */
  92.  
  93. struct edge_struct {
  94.     struct poly_struct *poly[2];    /* Each edge belongs to up to 2 polygons */
  95.     struct coordinate GPHUGE *vertex[2];    /* The two extreme points of this edge. */
  96.     struct edge_struct *next;    /* To chain lists */
  97.     int status,            /* Status flag to mark edges in scanning at certain Z level. */
  98.      position;            /* position in mesh: INNER_MESH, BOUNDARY or DIAGONNAL. */
  99. };
  100.  
  101. struct poly_struct {
  102.     struct edge_struct *edge[3];    /* As we do triangolation here... */
  103.     struct poly_struct *next;    /* To chain lists. */
  104. };
  105.  
  106. struct cntr_struct {        /* Contours are saved using this struct list. */
  107.     double X, Y;        /* The coordinates of this vertex. */
  108.     struct cntr_struct *next;    /* To chain lists. */
  109. };
  110.  
  111. static struct gnuplot_contours *contour_list = NULL;
  112. static double crnt_cntr[MAX_POINTS_PER_CNTR * 2];
  113. static int crnt_cntr_pt_index = 0;
  114. static double contour_level = 0.0;
  115. static int num_approx_pts = DEFAULT_NUM_APPROX_PTS;    /* # pts per approx/inter. */
  116. static int bspline_order = DEFAULT_BSPLINE_ORDER;    /* Bspline order to use. */
  117. static int interp_kind = INTERP_NOTHING;    /* Linear, Cubic interp., Bspline. */
  118. static double x_min, y_min, z_min;    /* Minimum values of x, y, and z */
  119. static double x_max, y_max, z_max;    /* Maximum values of x, y, and z */
  120.  
  121. static void add_cntr_point __PROTO((double x, double y));
  122. static void end_crnt_cntr __PROTO((void));
  123. static void gen_contours __PROTO((struct edge_struct * p_edges, double z_level,
  124.         double xx_min, double xx_max, double yy_min, double yy_max));
  125. static int update_all_edges __PROTO((struct edge_struct * p_edges,
  126.                      double z_level));
  127. static struct cntr_struct *gen_one_contour __PROTO((
  128.        struct edge_struct * p_edges, double z_level, int *contr_kind,
  129.                                int *num_active));
  130. static struct cntr_struct *trace_contour __PROTO((
  131.       struct edge_struct * pe_start, double z_level, int *num_active,
  132.                              int contr_kind));
  133. static struct cntr_struct *update_cntr_pt __PROTO((struct edge_struct * p_edge,
  134.                            double z_level));
  135. static int fuzzy_equal __PROTO((struct cntr_struct * p_cntr1,
  136.                 struct cntr_struct * p_cntr2));
  137.  
  138.  
  139. static void gen_triangle __PROTO((int num_isolines,
  140.          struct iso_curve * iso_lines, struct poly_struct ** p_polys,
  141.                   struct edge_struct ** p_edges));
  142. static void calc_min_max __PROTO((int num_isolines,
  143.                   struct iso_curve * iso_lines, double *xx_min, double *yy_min, double *zz_min,
  144.             double *xx_max, double *yy_max, double *zz_max));
  145. static struct edge_struct *add_edge __PROTO((struct coordinate GPHUGE * point0,
  146.      struct coordinate GPHUGE * point1, struct edge_struct ** p_edge,
  147.                      struct edge_struct ** pe_tail));
  148. static struct poly_struct *add_poly __PROTO((struct edge_struct * edge0,
  149.           struct edge_struct * edge1, struct edge_struct * edge2,
  150.        struct poly_struct ** p_poly, struct poly_struct ** pp_tail));
  151.  
  152.  
  153. static void put_contour __PROTO((struct cntr_struct * p_cntr, double z_level,
  154.           double xx_min, double xx_max, double yy_min, double yy_max,
  155.                  int contr_kind));
  156. static void put_contour_nothing __PROTO((struct cntr_struct * p_cntr));
  157. static int chk_contour_kind __PROTO((struct cntr_struct * p_cntr,
  158.                      int contr_kind));
  159. static void put_contour_cubic __PROTO((struct cntr_struct * p_cntr,
  160.                        double z_level, double xx_min, double xx_max, double yy_min, double yy_max,
  161.                        int contr_kind));
  162. static void put_contour_bspline __PROTO((struct cntr_struct * p_cntr,
  163.                      double z_level, double xx_min, double xx_max, double yy_min, double yy_max,
  164.                      int contr_kind));
  165. static void free_contour __PROTO((struct cntr_struct * p_cntr));
  166. static int count_contour __PROTO((struct cntr_struct * p_cntr));
  167. static int gen_cubic_spline __PROTO((int num_pts, struct cntr_struct * p_cntr,
  168.         double d2x[], double d2y[], double delta_t[], int contr_kind,
  169.                      double unit_x, double unit_y));
  170. static void intp_cubic_spline __PROTO((int n, struct cntr_struct * p_cntr,
  171.         double d2x[], double d2y[], double delta_t[], int n_intpol));
  172. static int solve_cubic_1 __PROTO((tri_diag m[], int n));
  173. static void solve_cubic_2 __PROTO((tri_diag m[], double x[], int n));
  174. /* 
  175.  * static int solve_tri_diag __PROTO((tri_diag m[], double r[], double x[], 
  176.  * int n));  see "protos.h"
  177.  */
  178. static void gen_bspline_approx __PROTO((struct cntr_struct * p_cntr,
  179.               int num_of_points, int order, int contr_kind));
  180. static void eval_bspline __PROTO((double t, struct cntr_struct * p_cntr,
  181.       int num_of_points, int order, int j, int contr_kind, double *x,
  182.                   double *y));
  183. static double fetch_knot __PROTO((int contr_kind, int num_of_points,
  184.                   int order, int i));
  185.  
  186. /*
  187.  * Entry routine to this whole set of contouring module.
  188.  */
  189. struct gnuplot_contours *contour(num_isolines, iso_lines, ZLevels, approx_pts, int_kind, order1, contour_levels_kind, cont_levels_list)
  190. int num_isolines;
  191. struct iso_curve *iso_lines;
  192. int ZLevels, approx_pts, int_kind, order1, contour_levels_kind;
  193. double *cont_levels_list;
  194. {
  195.     int i;
  196.     int num_of_z_levels;    /* # Z contour levels. */
  197.     struct poly_struct *p_polys, *p_poly;
  198.     struct edge_struct *p_edges, *p_edge;
  199.     double z = 0, dz = 0;
  200.     struct gnuplot_contours *save_contour_list;
  201.  
  202.     num_of_z_levels = ZLevels;
  203.     num_approx_pts = approx_pts;
  204.     bspline_order = order1 - 1;
  205.     interp_kind = int_kind;
  206.  
  207.     contour_list = NULL;
  208.  
  209.     /* 
  210.      * Calculate min/max values :
  211.      */
  212.     calc_min_max(num_isolines, iso_lines,
  213.          &x_min, &y_min, &z_min, &x_max, &y_max, &z_max);
  214.  
  215.     /* 
  216.      * Generate list of edges (p_edges) and list of triangles (p_polys):
  217.      */
  218.     gen_triangle(num_isolines, iso_lines, &p_polys, &p_edges);
  219.     crnt_cntr_pt_index = 0;
  220.  
  221.     if (contour_levels_kind == LEVELS_AUTO) {
  222.     dz = fabs(z_max - z_min);
  223.     if (dz == 0)
  224.         return NULL;    /* empty z range ? */
  225.     dz = set_tic(log10(dz), ((int) ZLevels + 1) * 2);
  226.     z = floor(z_min / dz) * dz;
  227.     num_of_z_levels = (int) floor((z_max - z) / dz);
  228.     }
  229.     for (i = 0; i < num_of_z_levels; i++) {
  230.     switch (contour_levels_kind) {
  231.     case LEVELS_AUTO:
  232.         z += dz;
  233.         break;
  234.     case LEVELS_INCREMENTAL:
  235.         z = cont_levels_list[0] + i * cont_levels_list[1];
  236.         break;
  237.     case LEVELS_DISCRETE:
  238.         z = is_log_z ? log(cont_levels_list[i]) / log_base_log_z : cont_levels_list[i];
  239.         break;
  240.     }
  241.     contour_level = z;
  242.     save_contour_list = contour_list;
  243.     gen_contours(p_edges, z, x_min, x_max, y_min, y_max);
  244.     if (contour_list != save_contour_list) {
  245.         contour_list->isNewLevel = 1;
  246.         sprintf(contour_list->label, contour_format, is_log_z ? pow(base_log_z, z) : z);
  247.     }
  248.     }
  249.  
  250.     /* Free all contouring related temporary data. */
  251.     while (p_polys) {
  252.     p_poly = p_polys->next;
  253.     free(p_polys);
  254.     p_polys = p_poly;
  255.     }
  256.     while (p_edges) {
  257.     p_edge = p_edges->next;
  258.     free(p_edges);
  259.     p_edges = p_edge;
  260.     }
  261.  
  262.     return contour_list;
  263. }
  264.  
  265. /*
  266.  * Adds another point to the currently build contour.
  267.  */
  268. static void add_cntr_point(x, y)
  269. double x, y;
  270. {
  271.     int index;
  272.  
  273.     if (crnt_cntr_pt_index >= MAX_POINTS_PER_CNTR - 1) {
  274.     index = crnt_cntr_pt_index - 1;
  275.     end_crnt_cntr();
  276.     crnt_cntr[0] = crnt_cntr[index * 2];
  277.     crnt_cntr[1] = crnt_cntr[index * 2 + 1];
  278.     crnt_cntr_pt_index = 1;    /* Keep the last point as first of this one. */
  279.     }
  280.     crnt_cntr[crnt_cntr_pt_index * 2] = x;
  281.     crnt_cntr[crnt_cntr_pt_index * 2 + 1] = y;
  282.     crnt_cntr_pt_index++;
  283. }
  284.  
  285. /*
  286.  * Done with current contour - create gnuplot data structure for it.
  287.  */
  288. static void end_crnt_cntr()
  289. {
  290.     int i;
  291.     struct gnuplot_contours *cntr = (struct gnuplot_contours *)
  292.     gp_alloc((unsigned long) sizeof(struct gnuplot_contours), "gnuplot_contour");
  293.     cntr->coords = (struct coordinate GPHUGE *)
  294.     gp_alloc((unsigned long) sizeof(struct coordinate)
  295.          * (unsigned long) crnt_cntr_pt_index, "contour coords");
  296.  
  297.     for (i = 0; i < crnt_cntr_pt_index; i++) {
  298.     cntr->coords[i].x = crnt_cntr[i * 2];
  299.     cntr->coords[i].y = crnt_cntr[i * 2 + 1];
  300.     cntr->coords[i].z = contour_level;
  301.     }
  302.     cntr->num_pts = crnt_cntr_pt_index;
  303.  
  304.     cntr->next = contour_list;
  305.     contour_list = cntr;
  306.     contour_list->isNewLevel = 0;
  307.  
  308.     crnt_cntr_pt_index = 0;
  309. }
  310.  
  311. /*
  312.  * Generates all contours by tracing the intersecting triangles.
  313.  */
  314. static void gen_contours(p_edges, z_level, xx_min, xx_max, yy_min, yy_max)
  315. struct edge_struct *p_edges;
  316. double z_level, xx_min, xx_max, yy_min, yy_max;
  317. {
  318.     int num_active,        /* Number of edges marked ACTIVE. */
  319.      contr_kind;        /* One of OPEN_CONTOUR, CLOSED_CONTOUR. */
  320.     struct cntr_struct *p_cntr;
  321.  
  322.     num_active = update_all_edges(p_edges, z_level);    /* Do pass 1. */
  323.  
  324.     contr_kind = OPEN_CONTOUR;    /* Start to look for contour on boundaries. */
  325.  
  326.     while (num_active > 0) {    /* Do Pass 2. */
  327.     /* Generate One contour (and update MumActive as needed): */
  328.     p_cntr = gen_one_contour(p_edges, z_level, &contr_kind, &num_active);
  329.     /* Emit it in requested format: */
  330.     put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind);
  331.     }
  332. }
  333.  
  334. /*
  335.  * Does pass 1, or marks the edges which are active (crosses this z_level)
  336.  * as ACTIVE, and the others as INACTIVE:
  337.  * Returns number of active edges (marked ACTIVE).
  338.  */
  339. static int update_all_edges(p_edges, z_level)
  340. struct edge_struct *p_edges;
  341. double z_level;
  342. {
  343.     int count = 0;
  344.  
  345.     while (p_edges) {
  346.     /* use the same test at both vertices to avoid roundoff errors */
  347.     if ((p_edges->vertex[0]->z >= z_level) !=
  348.         (p_edges->vertex[1]->z >= z_level)) {
  349.         p_edges->status = ACTIVE;
  350.         count++;
  351.     } else
  352.         p_edges->status = INACTIVE;
  353.     p_edges = p_edges->next;
  354.     }
  355.  
  356.     return count;
  357. }
  358.  
  359. /*
  360.  * Does pass 2, or find one complete contour out of the triangulation
  361.  * data base:
  362.  * Returns a pointer to the contour (as linked list), contr_kind is set to
  363.  * one of OPEN_CONTOUR, CLOSED_CONTOUR, and num_active is updated.
  364.  */
  365. static struct cntr_struct *gen_one_contour(p_edges, z_level, contr_kind, num_active)
  366. struct edge_struct *p_edges;    /* list of edges input */
  367. double z_level;    /* Z level of contour input */
  368. int *contr_kind;    /* OPEN_ or CLOESED_CONTOUR  in/out */
  369. int *num_active; /* number of active edges     in/out */
  370. {
  371.     struct edge_struct *pe_temp;
  372.  
  373.     if (*contr_kind == OPEN_CONTOUR) {
  374.     /* Look for something to start with on boundary: */
  375.     pe_temp = p_edges;
  376.     while (pe_temp) {
  377.         if ((pe_temp->status == ACTIVE) && (pe_temp->position == BOUNDARY))
  378.         break;
  379.         pe_temp = pe_temp->next;
  380.     }
  381.     if (!pe_temp)
  382.         *contr_kind = CLOSED_CONTOUR;    /* No more contours on boundary. */
  383.     else {
  384.         return trace_contour(pe_temp, z_level, num_active, *contr_kind);
  385.     }
  386.     }
  387.     if (*contr_kind == CLOSED_CONTOUR) {
  388.     /* Look for something to start with inside: */
  389.     pe_temp = p_edges;
  390.     while (pe_temp) {
  391.         if ((pe_temp->status == ACTIVE) && (!(pe_temp->position == BOUNDARY)))
  392.         break;
  393.         pe_temp = pe_temp->next;
  394.     }
  395.     if (!pe_temp) {
  396.         *num_active = 0;
  397.         fprintf(stderr, "gen_one_contour: no contour found\n");
  398.         return NULL;
  399.     } else {
  400.         *contr_kind = CLOSED_CONTOUR;
  401.         return trace_contour(pe_temp, z_level, num_active, *contr_kind);
  402.     }
  403.     }
  404.     return NULL;        /* We should never be here, but lint... */
  405. }
  406.  
  407. /*
  408.  * Search the data base along a contour starts at the edge pe_start until
  409.  * a boundary edge is detected or until we close the loop back to pe_start.
  410.  * Returns a linked list of all the points on the contour
  411.  * Also decreases num_active by the number of points on contour.
  412.  */
  413. static struct cntr_struct *trace_contour(pe_start, z_level, num_active, contr_kind)
  414. struct edge_struct *pe_start;    /* edge to start contour input */
  415. double z_level;    /* Z level of contour input */
  416. int *num_active;    /* number of active edges in/out */
  417. int contr_kind;        /* OPEN_ or CLOESED_CONTOUR    input */
  418. {
  419.     struct cntr_struct *p_cntr, *pc_tail;
  420.     struct edge_struct *p_edge, *p_next_edge;
  421.     struct poly_struct *p_poly, *PLastpoly = NULL;
  422.     int i;
  423.  
  424.     p_edge = pe_start;        /* first edge to start contour */
  425.  
  426.     /* Generate the header of the contour - the point on pe_start. */
  427.     if (contr_kind == OPEN_CONTOUR) {
  428.     pe_start->status = INACTIVE;
  429.     (*num_active)--;
  430.     }
  431.     if (p_edge->poly[0] || p_edge->poly[1]) {    /* more than one point */
  432.  
  433.     p_cntr = pc_tail = update_cntr_pt(pe_start, z_level);    /* first point */
  434.  
  435.     do {
  436.         /* Find polygon to continue (Not where we came from - PLastpoly): */
  437.         if (p_edge->poly[0] == PLastpoly)
  438.         p_poly = p_edge->poly[1];
  439.         else
  440.         p_poly = p_edge->poly[0];
  441.         p_next_edge = NULL;    /* In case of error, remains NULL. */
  442.         for (i = 0; i < 3; i++)    /* Test the 3 edges of the polygon: */
  443.         if (p_poly->edge[i] != p_edge)
  444.             if (p_poly->edge[i]->status == ACTIVE)
  445.             p_next_edge = p_poly->edge[i];
  446.         if (!p_next_edge) {    /* Error exit */
  447.         pc_tail->next = NULL;
  448.         free_contour(p_cntr);
  449.         fprintf(stderr, "trace_contour: unexpected end of contour\n");
  450.         return NULL;
  451.         }
  452.         p_edge = p_next_edge;
  453.         PLastpoly = p_poly;
  454.         p_edge->status = INACTIVE;
  455.         (*num_active)--;
  456.  
  457.         /* Do not allocate contour points on diagonal edges */
  458.         if (p_edge->position != DIAGONAL) {
  459.  
  460.         pc_tail->next = update_cntr_pt(p_edge, z_level);
  461.  
  462.         /* Remove nearby points */
  463.         if (fuzzy_equal(pc_tail, pc_tail->next)) {
  464.  
  465.             free((char *) pc_tail->next);
  466.         } else
  467.             pc_tail = pc_tail->next;
  468.         }
  469.     } while ((p_edge != pe_start) && (p_edge->position != BOUNDARY));
  470.  
  471.     pc_tail->next = NULL;
  472.  
  473.     /* For CLOSED_CONTOUR the first and last point should be equal */
  474.     if (pe_start == p_edge) {
  475.         (p_cntr->X) = (pc_tail->X);
  476.         (p_cntr->Y) = (pc_tail->Y);
  477.     }
  478.     } else {            /* only one point, forget it */
  479.     p_cntr = NULL;
  480.     }
  481.  
  482.     return p_cntr;
  483. }
  484.  
  485. /*
  486.  * Allocates one contour location and update it to to correct position
  487.  * according to z_level and edge p_edge.
  488.  */
  489. static struct cntr_struct *update_cntr_pt(p_edge, z_level)
  490. struct edge_struct *p_edge;
  491. double z_level;
  492. {
  493.     double t;
  494.     struct cntr_struct *p_cntr;
  495.  
  496.     t = (z_level - p_edge->vertex[0]->z) /
  497.     (p_edge->vertex[1]->z - p_edge->vertex[0]->z);
  498.  
  499.     /* test if t is out of interval [0:1] (should not happen but who knows ...) */
  500.     t = (t < 0.0 ? 0.0 : t);
  501.     t = (t > 1.0 ? 1.0 : t);
  502.  
  503.     p_cntr = (struct cntr_struct *)
  504.     gp_alloc((unsigned long) sizeof(struct cntr_struct), "contour cntr_struct");
  505.  
  506.     p_cntr->X = p_edge->vertex[1]->x * t +
  507.     p_edge->vertex[0]->x * (1 - t);
  508.     p_cntr->Y = p_edge->vertex[1]->y * t +
  509.     p_edge->vertex[0]->y * (1 - t);
  510.     return p_cntr;
  511. }
  512.  
  513. /*
  514.  * Simple routine to decide if two contour points are equal by
  515.  * calculating the relative error (< EPSILON).
  516.  */
  517. static int fuzzy_equal(p_cntr1, p_cntr2)
  518. struct cntr_struct *p_cntr1, *p_cntr2;
  519. {
  520.     double unit_x, unit_y;
  521.     unit_x = ABS(x_max - x_min) + zero;        /* reference */
  522.     unit_y = ABS(y_max - y_min) + zero;
  523.     return (
  524.            ABS(p_cntr1->X - p_cntr2->X) / unit_x < EPSILON &&
  525.            ABS(p_cntr1->Y - p_cntr2->Y) / unit_y < EPSILON);
  526. }
  527.  
  528. /*
  529.  * Generate the triangles.
  530.  * Returns the lists (edges & polys) via pointers to their heads.
  531.  */
  532. static void gen_triangle(num_isolines, iso_lines, p_polys, p_edges)
  533. int num_isolines;    /* number of iso-lines input */
  534. struct iso_curve *iso_lines;    /* iso-lines input */
  535. struct poly_struct **p_polys;    /* list of polygons output */
  536. struct edge_struct **p_edges; /* list of edges output */
  537. {
  538.     int i, j, grid_x_max = iso_lines->p_count;
  539.     struct edge_struct *p_edge1, *p_edge2, *edge0, *edge1, *edge2, *pe_tail,
  540.     *pe_tail1, *pe_tail2, *pe_temp;
  541.     struct poly_struct *pp_tail, *lower_tri, *upper_tri;
  542.     struct coordinate GPHUGE *p_vrtx1, GPHUGE * p_vrtx2;    /* HBB 980308: need to tag *each* of them as GPHUGE! */
  543.  
  544.     (*p_polys) = pp_tail = NULL;    /* clear lists */
  545.     (*p_edges) = pe_tail = NULL;
  546.  
  547.     p_vrtx1 = iso_lines->points;    /* first row of vertices */
  548.     p_edge1 = pe_tail1 = NULL;    /* clear list of edges */
  549.  
  550.     /* Generate edges of first row */
  551.     for (j = 0; j < grid_x_max - 1; j++)
  552.     add_edge(p_vrtx1 + j, p_vrtx1 + j + 1, &p_edge1, &pe_tail1);
  553.  
  554.     (*p_edges) = p_edge1;    /* update main list */
  555.     pe_tail = pe_tail1;
  556.  
  557.  
  558.     /*
  559.      * Combines vertices to edges and edges to triangles:
  560.      * ==================================================
  561.      * The edges are stored in the edge list, referenced by p_edges
  562.      * (pe_tail points on last edge).
  563.      *
  564.      * Temporary pointers:
  565.      * 1. p_edge2: Top horizontal edge list:       -----------------------  2
  566.      * 2. pe_tail: middle edge list:              |\  |\  |\  |\  |\  |\  |
  567.      *                                            |  \|  \|  \|  \|  \|  \|
  568.      * 3. p_edge1: Bottom horizontal edge list:    -----------------------  1
  569.      *
  570.      * The routine generates two triangle            Lower      Upper 1  
  571.      * upper one and lower one:                     | \           ----   
  572.      * (Nums. are edges order in polys)            0|   \1       0\   |2 
  573.      * The polygons are stored in the polygon        ----           \ |  
  574.      * list (*p_polys) (pp_tail points on             2                  
  575.      * last polygon).
  576.      *                                                        1
  577.      *                                                   -----------
  578.      * In addition, the edge lists are updated -        | \   0     |
  579.      * each edge has two pointers on the two            |   \       |
  580.      * (one active if boundary) polygons which         0|1   0\1   0|1
  581.      * uses it. These two pointer to polygons           |       \   |
  582.      * are named: poly[0], poly[1]. The diagram         |    1    \ |
  583.      * on the right show how they are used for the       -----------
  584.      * upper and lower polygons (INNER_MESH polygons only).  0
  585.      */
  586.  
  587.     for (i = 1; i < num_isolines; i++) {
  588.     /* Read next column and gen. polys. */
  589.     iso_lines = iso_lines->next;
  590.  
  591.     p_vrtx2 = iso_lines->points;    /* next row of vertices */
  592.     p_edge2 = pe_tail2 = NULL;    /* clear top horizontal list */
  593.     pe_temp = p_edge1;    /* pointer in bottom list */
  594.  
  595.     /*
  596.      * Generate edges and triagles for next row:
  597.      */
  598.  
  599.     /* generate first vertical edge */
  600.     edge2 = add_edge(p_vrtx1, p_vrtx2, p_edges, &pe_tail);
  601.  
  602.     for (j = 0; j < grid_x_max - 1; j++) {
  603.  
  604.         /* copy vertical edge for lower triangle */
  605.         edge0 = edge2;
  606.  
  607.         if (pe_temp && pe_temp->vertex[0] == p_vrtx1 + j) {
  608.         /* test lower edge */
  609.         edge2 = pe_temp;
  610.         pe_temp = pe_temp->next;
  611.         } else {
  612.         edge2 = NULL;    /* edge is undefined */
  613.         }
  614.  
  615.         /* generate diagonal edge */
  616.         edge1 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j, p_edges, &pe_tail);
  617.         if (edge1)
  618.         edge1->position = DIAGONAL;
  619.  
  620.         /* generate lower triangle */
  621.         lower_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
  622.  
  623.         /* copy diagonal edge for upper triangle */
  624.         edge0 = edge1;
  625.  
  626.         /* generate upper edge */
  627.         edge1 = add_edge(p_vrtx2 + j, p_vrtx2 + j + 1, &p_edge2, &pe_tail2);
  628.  
  629.         /* generate vertical edge */
  630.         edge2 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j + 1, p_edges, &pe_tail);
  631.  
  632.         /* generate upper triangle */
  633.         upper_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
  634.     }
  635.  
  636.     if ((*p_edges)) {    /* Chain new edges to main list. */
  637.         pe_tail->next = p_edge2;
  638.         pe_tail = pe_tail2;
  639.     } else {
  640.         (*p_edges) = p_edge2;
  641.         pe_tail = pe_tail2;
  642.     }
  643.  
  644.     p_edge1 = p_edge2;
  645.     p_vrtx1 = p_vrtx2;
  646.     }
  647.  
  648.     /* Update the boundary flag, saved in each edge, and update indexes: */
  649.  
  650.     pe_temp = (*p_edges);
  651.  
  652.     while (pe_temp) {
  653.     if ((!(pe_temp->poly[0])) || (!(pe_temp->poly[1])))
  654.         (pe_temp->position) = BOUNDARY;
  655.     pe_temp = pe_temp->next;
  656.     }
  657. }
  658.  
  659. /*
  660.  * Calculate minimum and maximum values
  661.  */
  662. static void calc_min_max(num_isolines, iso_lines, xx_min, yy_min, zz_min, xx_max, yy_max, zz_max)
  663. int num_isolines;    /* number of iso-lines input */
  664. struct iso_curve *iso_lines;    /* iso-lines input */
  665. double *xx_min, *yy_min, *zz_min, *xx_max, *yy_max, *zz_max;  /* min/max values in/out */
  666. {
  667.     int i, j, grid_x_max;
  668.     struct coordinate GPHUGE *vertex;
  669.  
  670.     grid_x_max = iso_lines->p_count;    /* number of vertices per iso_line */
  671.  
  672.     (*xx_min) = (*yy_min) = (*zz_min) = VERYLARGE; /* clear min/max values */
  673.     (*xx_max) = (*yy_max) = (*zz_max) = -VERYLARGE;
  674.  
  675.     for (j = 0; j < num_isolines; j++) {
  676.  
  677.     vertex = iso_lines->points;
  678.  
  679.     for (i = 0; i < grid_x_max; i++) {
  680.         if (vertex[i].type != UNDEFINED) {
  681.         if (vertex[i].x > (*xx_max))
  682.             (*xx_max) = vertex[i].x;
  683.         if (vertex[i].y > (*yy_max))
  684.             (*yy_max) = vertex[i].y;
  685.         if (vertex[i].z > (*zz_max))
  686.             (*zz_max) = vertex[i].z;
  687.         if (vertex[i].x < (*xx_min))
  688.             (*xx_min) = vertex[i].x;
  689.         if (vertex[i].y < (*yy_min))
  690.             (*yy_min) = vertex[i].y;
  691.         if (vertex[i].z < (*zz_min))
  692.             (*zz_min) = vertex[i].z;
  693.         }
  694.     }
  695.     iso_lines = iso_lines->next;
  696.     }
  697.     /* 
  698.      * fprintf(stderr," x: %g, %g\n", (*xx_min), (*xx_max));
  699.      * fprintf(stderr," y: %g, %g\n", (*yy_min), (*yy_max));
  700.      * fprintf(stderr," z: %g, %g\n", (*zz_min), (*zz_max));
  701.      */
  702. }
  703.  
  704. /*
  705.  * Generate new edge and append it to list, but only if both vertices are 
  706.  * defined. The list is referenced by p_edge and pe_tail (p_edge points on 
  707.  * first edge and pe_tail on last one).
  708.  * Note, the list may be empty (pe_edge==pe_tail==NULL) on entry and exit.
  709.  */
  710. static struct edge_struct *add_edge(point0, point1, p_edge, pe_tail)
  711. struct coordinate GPHUGE * point0;    /* 2 vertices input */
  712. struct coordinate GPHUGE * point1;
  713. struct edge_struct **p_edge, **pe_tail;  /* pointers to edge list in/out */
  714. {
  715.     struct edge_struct *pe_temp = NULL;
  716.  
  717.     if (point0->type != UNDEFINED && point1->type != UNDEFINED) {
  718.  
  719.     pe_temp = (struct edge_struct *)
  720.         gp_alloc((unsigned long) sizeof(struct edge_struct), "contour edge");
  721.  
  722.     pe_temp->poly[0] = NULL;    /* clear links           */
  723.     pe_temp->poly[1] = NULL;
  724.     pe_temp->vertex[0] = point0;    /* First vertex of edge. */
  725.     pe_temp->vertex[1] = point1;    /* Second vertex of edge. */
  726.     pe_temp->next = NULL;
  727.     pe_temp->position = INNER_MESH;        /* default position in mesh */
  728.  
  729.     if ((*pe_tail)) {
  730.         (*pe_tail)->next = pe_temp;        /* Stick new record as last one. */
  731.     } else {
  732.         (*p_edge) = pe_temp;    /* start new list if empty */
  733.     }
  734.     (*pe_tail) = pe_temp;    /* continue to last record. */
  735.  
  736.     }
  737.     return pe_temp;        /* returns NULL, if no edge allocated */
  738. }
  739.  
  740. /*
  741.  * Generate new triangle and append it to list, but only if all edges are defined.
  742.  * The list is referenced by p_poly and pp_tail (p_poly points on first ploygon
  743.  * and pp_tail on last one).
  744.  * Note, the list may be empty (pe_ploy==pp_tail==NULL) on entry and exit.
  745.  */
  746. static struct poly_struct *add_poly(edge0, edge1, edge2, p_poly, pp_tail)
  747. struct edge_struct *edge0, *edge1, *edge2;    /* 3 edges input */
  748. struct poly_struct **p_poly, **pp_tail;    /* pointers to polygon list in/out */
  749. {
  750.     struct poly_struct *pp_temp = NULL;
  751.  
  752.     if (edge0 && edge1 && edge2) {
  753.  
  754.     pp_temp = (struct poly_struct *)
  755.         gp_alloc((unsigned long) sizeof(struct poly_struct), "contour polygon");
  756.  
  757.     pp_temp->edge[0] = edge0;    /* First edge of triangle */
  758.     pp_temp->edge[1] = edge1;    /* Second one             */
  759.     pp_temp->edge[2] = edge2;    /* Third one              */
  760.     pp_temp->next = NULL;
  761.  
  762.     if (edge0->poly[0])    /* update edge0 */
  763.         edge0->poly[1] = pp_temp;
  764.     else
  765.         edge0->poly[0] = pp_temp;
  766.  
  767.     if (edge1->poly[0])    /* update edge1 */
  768.         edge1->poly[1] = pp_temp;
  769.     else
  770.         edge1->poly[0] = pp_temp;
  771.  
  772.     if (edge2->poly[0])    /* update edge2 */
  773.         edge2->poly[1] = pp_temp;
  774.     else
  775.         edge2->poly[0] = pp_temp;
  776.  
  777.     if ((*pp_tail))        /* Stick new record as last one. */
  778.         (*pp_tail)->next = pp_temp;
  779.     else
  780.         (*p_poly) = pp_temp;    /* start new list if empty */
  781.  
  782.     (*pp_tail) = pp_temp;    /* continue to last record. */
  783.  
  784.     }
  785.     return pp_temp;        /* returns NULL, if no edge allocated */
  786. }
  787.  
  788.  
  789.  
  790. /*
  791.  * Calls the (hopefully) desired interpolation/approximation routine.
  792.  */
  793. static void put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind)
  794. struct cntr_struct *p_cntr;    /* contour structure input */
  795. double z_level,    /* Z level of contour input */
  796.   xx_min, xx_max, yy_min, yy_max; /* minimum/maximum values input */
  797. int contr_kind;            /* OPEN_ or CLOESED_CONTOUR input */
  798. {
  799.  
  800.     if (!p_cntr)
  801.     return;            /* Nothing to do if it is empty contour. */
  802.  
  803.     switch (interp_kind) {
  804.     case INTERP_NOTHING:    /* No interpolation/approximation. */
  805.     put_contour_nothing(p_cntr);
  806.     break;
  807.     case INTERP_CUBIC:        /* Cubic spline interpolation. */
  808.     put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
  809.               chk_contour_kind(p_cntr, contr_kind));
  810.  
  811.     break;
  812.     case APPROX_BSPLINE:    /* Bspline approximation. */
  813.     put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
  814.                 chk_contour_kind(p_cntr, contr_kind));
  815.     break;
  816.     }
  817.     free_contour(p_cntr);
  818. }
  819.  
  820. /*
  821.  * Simply puts contour coordinates in order with no interpolation or
  822.  * approximation.
  823.  */
  824. static void put_contour_nothing(p_cntr)
  825. struct cntr_struct *p_cntr;
  826. {
  827.     while (p_cntr) {
  828.     add_cntr_point(p_cntr->X, p_cntr->Y);
  829.     p_cntr = p_cntr->next;
  830.     }
  831.     end_crnt_cntr();
  832. }
  833.  
  834. /*
  835.  * for some reason contours are never flagged as CLOSED_CONTOUR
  836.  * if first point == last point, set flag accordingly
  837.  *
  838.  */
  839.  
  840. static int chk_contour_kind(p_cntr, contr_kind)
  841. struct cntr_struct *p_cntr;
  842. int contr_kind;
  843. {
  844.     struct cntr_struct *pc_tail = NULL;
  845.     int current_contr_kind;
  846.  
  847.     FPRINTF((stderr, "check_contour_kind: current contr_kind value is %d\n", contr_kind));
  848.  
  849.     current_contr_kind = contr_kind;
  850.  
  851.     if (contr_kind != CLOSED_CONTOUR) {
  852.     pc_tail = p_cntr;
  853.     while (pc_tail->next)
  854.         pc_tail = pc_tail->next;    /* Find last point. */
  855.  
  856.     /* test if first and last point are equal */
  857.     if (fuzzy_equal(pc_tail, p_cntr)) {
  858.         current_contr_kind = CLOSED_CONTOUR;
  859.         FPRINTF((stderr, "check_contour_kind: contr_kind changed to %d\n", current_contr_kind));
  860.     }
  861.     }
  862.     return (current_contr_kind);
  863. }
  864.  
  865. /*
  866.  * Generate a cubic spline curve through the points (x_i,y_i) which are
  867.  * stored in the linked list p_cntr.
  868.  * The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
  869.  * parameter t is the length of the linear stroke.
  870.  */
  871. static void put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind)
  872. struct cntr_struct *p_cntr;
  873. double z_level, xx_min, xx_max, yy_min, yy_max;
  874. int contr_kind;
  875. {
  876.     int num_pts, num_intpol;
  877.     double unit_x, unit_y;    /* To define norm (x,y)-plane */
  878.     double *delta_t;        /* Interval length t_{i+1}-t_i */
  879.     double *d2x, *d2y;        /* Second derivatives x''(t_i), y''(t_i) */
  880.     struct cntr_struct *pc_tail;
  881.  
  882.     num_pts = count_contour(p_cntr);    /* Number of points in contour. */
  883.  
  884.     pc_tail = p_cntr;        /* Find last point. */
  885.     while (pc_tail->next)
  886.     pc_tail = pc_tail->next;
  887.  
  888.     if (contr_kind == CLOSED_CONTOUR) {
  889.     /* Test if first and last point are equal (should be) */
  890.     if (!fuzzy_equal(pc_tail, p_cntr)) {
  891.         pc_tail->next = p_cntr;    /* Close contour list - make it circular. */
  892.         num_pts++;
  893.     }
  894.     }
  895.     delta_t = (double *)
  896.     gp_alloc((unsigned long) (sizeof(double) * num_pts), "contour delta_t");
  897.     d2x = (double *)
  898.     gp_alloc((unsigned long) (sizeof(double) * num_pts), "contour d2x");
  899.     d2y = (double *)
  900.     gp_alloc((unsigned long) (sizeof(double) * num_pts), "contour d2y");
  901.  
  902.     /* Width and hight of the grid is used at unit length (2d-norm) */
  903.     unit_x = xx_max - x_min;
  904.     unit_y = yy_max - y_min;
  905.     unit_x = (unit_x > zero ? unit_x : zero);    /* should not be zero */
  906.     unit_y = (unit_y > zero ? unit_y : zero);
  907.  
  908.     if (num_pts > 2) {
  909.     /*
  910.      * Calculate second derivatives d2x[], d2y[] and interval lengths delta_t[]:
  911.      */
  912.     if (!gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t,
  913.                   contr_kind, unit_x, unit_y)) {
  914.         free((char *) delta_t);
  915.         free((char *) d2x);
  916.         free((char *) d2y);
  917.         if (contr_kind == CLOSED_CONTOUR)
  918.         pc_tail->next = NULL;    /* Un-circular list */
  919.         return;
  920.     }
  921.     }
  922.     /* If following (num_pts > 1) is TRUE then exactly 2 points in contour.  */
  923.     else if (num_pts > 1) {
  924.     /* set all second derivatives to zero, interval length to 1 */
  925.     d2x[0] = 0.;
  926.     d2y[0] = 0.;
  927.     d2x[1] = 0.;
  928.     d2y[1] = 0.;
  929.     delta_t[0] = 1.;
  930.     } else {            /* Only one point ( ?? ) - ignore it. */
  931.     free((char *) delta_t);
  932.     free((char *) d2x);
  933.     free((char *) d2y);
  934.     if (contr_kind == CLOSED_CONTOUR)
  935.         pc_tail->next = NULL;    /* Un-circular list */
  936.     return;
  937.     }
  938.  
  939.     /* Calculate "num_intpol" interpolated values */
  940.     num_intpol = 1 + (num_pts - 1) * num_approx_pts;    /* global: num_approx_pts */
  941.     intp_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, num_intpol);
  942.  
  943.     free((char *) delta_t);
  944.     free((char *) d2x);
  945.     free((char *) d2y);
  946.  
  947.     if (contr_kind == CLOSED_CONTOUR)
  948.     pc_tail->next = NULL;    /* Un-circular list */
  949.  
  950.     end_crnt_cntr();
  951. }
  952.  
  953.  
  954. /*
  955.  * Find Bspline approximation for this data set.
  956.  * Uses global variable num_approx_pts to determine number of samples per
  957.  * interval, where the knot vector intervals are assumed to be uniform, and
  958.  * Global variable bspline_order for the order of Bspline to use.
  959.  */
  960. static void put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind)
  961. struct cntr_struct *p_cntr;
  962. double z_level, xx_min, xx_max, yy_min, yy_max;
  963. int contr_kind;
  964. {
  965.     int num_pts, order = bspline_order;
  966.  
  967.     num_pts = count_contour(p_cntr);    /* Number of points in contour. */
  968.     if (num_pts < 2)
  969.     return;            /* Can't do nothing if empty or one points! */
  970.     /* Order must be less than number of points in curve - fix it if needed. */
  971.     if (order > num_pts - 1)
  972.     order = num_pts - 1;
  973.  
  974.     gen_bspline_approx(p_cntr, num_pts, order, contr_kind);
  975.     end_crnt_cntr();
  976. }
  977.  
  978. /*
  979.  * Free all elements in the contour list.
  980.  */
  981. static void free_contour(p_cntr)
  982. struct cntr_struct *p_cntr;
  983. {
  984.     struct cntr_struct *pc_temp;
  985.  
  986.     while (p_cntr) {
  987.     pc_temp = p_cntr;
  988.     p_cntr = p_cntr->next;
  989.     free((char *) pc_temp);
  990.     }
  991. }
  992.  
  993. /*
  994.  * Counts number of points in contour.
  995.  */
  996. static int count_contour(p_cntr)
  997. struct cntr_struct *p_cntr;
  998. {
  999.     int count = 0;
  1000.  
  1001.     while (p_cntr) {
  1002.     count++;
  1003.     p_cntr = p_cntr->next;
  1004.     }
  1005.     return count;
  1006. }
  1007.  
  1008. /*
  1009.  * Find second derivatives (x''(t_i),y''(t_i)) of cubic spline interpolation
  1010.  * through list of points (x_i,y_i). The parameter t is calculated as the
  1011.  * length of the linear stroke. The number of points must be at least 3.
  1012.  * Note: For CLOSED_CONTOURs the first and last point must be equal.
  1013.  */
  1014. static int gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, contr_kind, unit_x, unit_y)
  1015. int num_pts;    /* Number of points (num_pts>=3), input */
  1016. struct cntr_struct *p_cntr;    /* List of points (x(t_i),y(t_i)), input */
  1017. double d2x[], d2y[],    /* Second derivatives (x''(t_i),y''(t_i)), output */
  1018.   delta_t[];    /* List of interval lengths t_{i+1}-t_{i}, output */
  1019. int contr_kind;    /* CLOSED_CONTOUR or OPEN_CONTOUR, input  */
  1020. double unit_x, unit_y;    /* Unit length in x and y (norm=1), input */
  1021. {
  1022.     int n, i;
  1023.     double norm;
  1024.     tri_diag *m;        /* The tri-diagonal matrix is saved here. */
  1025.     struct cntr_struct *pc_temp;
  1026.  
  1027.     m = (tri_diag *)
  1028.     gp_alloc((unsigned long) (sizeof(tri_diag) * num_pts), "contour tridiag m");
  1029.  
  1030.     /*
  1031.      * Calculate first differences in (d2x[i], d2y[i]) and interval lengths
  1032.      * in delta_t[i]:
  1033.      */
  1034.     pc_temp = p_cntr;
  1035.     for (i = 0; i < num_pts - 1; i++) {
  1036.     d2x[i] = pc_temp->next->X - pc_temp->X;
  1037.     d2y[i] = pc_temp->next->Y - pc_temp->Y;
  1038.     /*
  1039.      * The Norm of a linear stroke is calculated in "normal coordinates"
  1040.      * and used as interval length:
  1041.      */
  1042.     delta_t[i] = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y));
  1043.  
  1044.     d2x[i] /= delta_t[i];    /* first difference, with unit norm: */
  1045.     d2y[i] /= delta_t[i];    /*   || (d2x[i], d2y[i]) || = 1      */
  1046.  
  1047.     pc_temp = pc_temp->next;
  1048.     }
  1049.  
  1050.     /*
  1051.      * Setup linear System:  M * x = b
  1052.      */
  1053.     n = num_pts - 2;        /* Without first and last point */
  1054.     if (contr_kind == CLOSED_CONTOUR) {
  1055.     /* First and last points must be equal for CLOSED_CONTOURs */
  1056.     delta_t[num_pts - 1] = delta_t[0];
  1057.     d2x[num_pts - 1] = d2x[0];
  1058.     d2y[num_pts - 1] = d2y[0];
  1059.     n++;            /* Add last point (= first point) */
  1060.     }
  1061.     for (i = 0; i < n; i++) {
  1062.     /* Matrix M, mainly tridiagonal with cyclic second index ("j = j+n mod n") */
  1063.     m[i][0] = delta_t[i];    /* Off-diagonal element M_{i,i-1} */
  1064.     m[i][1] = 2. * (delta_t[i] + delta_t[i + 1]);    /* M_{i,i} */
  1065.     m[i][2] = delta_t[i + 1];    /* Off-diagonal element M_{i,i+1} */
  1066.  
  1067.     /* Right side b_x and b_y */
  1068.     d2x[i] = (d2x[i + 1] - d2x[i]) * 6.;
  1069.     d2y[i] = (d2y[i + 1] - d2y[i]) * 6.;
  1070.  
  1071.     /*
  1072.      * If the linear stroke shows a cusps of more than 90 degree, the right
  1073.      * side is reduced to avoid oscillations in the spline:
  1074.      */
  1075.     norm = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y)) / 8.5;
  1076.  
  1077.     if (norm > 1.) {
  1078.         d2x[i] /= norm;
  1079.         d2y[i] /= norm;
  1080.         /* The first derivative will not be continuous */
  1081.     }
  1082.     }
  1083.  
  1084.     if (contr_kind != CLOSED_CONTOUR) {
  1085.     /* Third derivative is set to zero at both ends */
  1086.     m[0][1] += m[0][0];    /* M_{0,0}     */
  1087.     m[0][0] = 0.;        /* M_{0,n-1}   */
  1088.     m[n - 1][1] += m[n - 1][2];    /* M_{n-1,n-1} */
  1089.     m[n - 1][2] = 0.;    /* M_{n-1,0}   */
  1090.     }
  1091.     /* Solve linear systems for d2x[] and d2y[] */
  1092.  
  1093.  
  1094.     if (solve_cubic_1(m, n)) {    /* Calculate Cholesky decomposition */
  1095.     solve_cubic_2(m, d2x, n);    /* solve M * d2x = b_x */
  1096.     solve_cubic_2(m, d2y, n);    /* solve M * d2y = b_y */
  1097.  
  1098.     } else {            /* Should not happen, but who knows ... */
  1099.     free((char *) m);
  1100.     return FALSE;
  1101.     }
  1102.  
  1103.     /* Shift all second derivatives one place right and abdate end points */
  1104.     for (i = n; i > 0; i--) {
  1105.     d2x[i] = d2x[i - 1];
  1106.     d2y[i] = d2y[i - 1];
  1107.     }
  1108.     if (contr_kind == CLOSED_CONTOUR) {
  1109.     d2x[0] = d2x[n];
  1110.     d2y[0] = d2y[n];
  1111.     } else {
  1112.     d2x[0] = d2x[1];    /* Third derivative is zero in */
  1113.     d2y[0] = d2y[1];    /*     first and last interval */
  1114.     d2x[n + 1] = d2x[n];
  1115.     d2y[n + 1] = d2y[n];
  1116.     }
  1117.  
  1118.     free((char *) m);
  1119.     return TRUE;
  1120. }
  1121.  
  1122. /*
  1123.  * Calculate interpolated values of the spline function (defined via p_cntr
  1124.  * and the second derivatives d2x[] and d2y[]). The number of tabulated
  1125.  * values is n. On an equidistant grid n_intpol values are calculated.
  1126.  */
  1127. static void intp_cubic_spline(n, p_cntr, d2x, d2y, delta_t, n_intpol)
  1128. int n;
  1129. struct cntr_struct *p_cntr;
  1130. double d2x[], d2y[], delta_t[];
  1131. int n_intpol;
  1132. {
  1133.     double t, t_skip, t_max;
  1134.     double x0, x1, x, y0, y1, y;
  1135.     double d, hx, dx0, dx01, hy, dy0, dy01;
  1136.     int i;
  1137.  
  1138.     /* The length of the total interval */
  1139.     t_max = 0.;
  1140.     for (i = 0; i < n - 1; i++)
  1141.     t_max += delta_t[i];
  1142.  
  1143.     /* The distance between interpolated points */
  1144.     t_skip = (1. - 1e-7) * t_max / (n_intpol - 1);
  1145.  
  1146.     t = 0.;            /* Parameter value */
  1147.     x1 = p_cntr->X;
  1148.     y1 = p_cntr->Y;
  1149.     add_cntr_point(x1, y1);    /* First point. */
  1150.     t += t_skip;
  1151.  
  1152.     for (i = 0; i < n - 1; i++) {
  1153.     p_cntr = p_cntr->next;
  1154.  
  1155.     d = delta_t[i];        /* Interval length */
  1156.     x0 = x1;
  1157.     y0 = y1;
  1158.     x1 = p_cntr->X;
  1159.     y1 = p_cntr->Y;
  1160.     hx = (x1 - x0) / d;
  1161.     hy = (y1 - y0) / d;
  1162.     dx0 = (d2x[i + 1] + 2 * d2x[i]) / 6.;
  1163.     dy0 = (d2y[i + 1] + 2 * d2y[i]) / 6.;
  1164.     dx01 = (d2x[i + 1] - d2x[i]) / (6. * d);
  1165.     dy01 = (d2y[i + 1] - d2y[i]) / (6. * d);
  1166.     while (t <= delta_t[i]) {    /* t in current interval ? */
  1167.         x = x0 + t * (hx + (t - d) * (dx0 + t * dx01));
  1168.         y = y0 + t * (hy + (t - d) * (dy0 + t * dy01));
  1169.         add_cntr_point(x, y);    /* next point. */
  1170.         t += t_skip;
  1171.     }
  1172.     t -= delta_t[i];    /* Parameter t relative to start of next interval */
  1173.     }
  1174. }
  1175.  
  1176. /*
  1177.  * The following two procedures solve the special linear system which arise
  1178.  * in cubic spline interpolation. If x is assumed cyclic ( x[i]=x[n+i] ) the
  1179.  * equations can be written as (i=0,1,...,n-1):
  1180.  *     m[i][0] * x[i-1] + m[i][1] * x[i] + m[i][2] * x[i+1] = b[i] .
  1181.  * In matrix notation one gets M * x = b, where the matrix M is tridiagonal
  1182.  * with additional elements in the upper right and lower left position:
  1183.  *   m[i][0] = M_{i,i-1}  for i=1,2,...,n-1    and    m[0][0] = M_{0,n-1} ,
  1184.  *   m[i][1] = M_{i, i }  for i=0,1,...,n-1
  1185.  *   m[i][2] = M_{i,i+1}  for i=0,1,...,n-2    and    m[n-1][2] = M_{n-1,0}.
  1186.  * M should be symmetric (m[i+1][0]=m[i][2]) and positiv definite.
  1187.  * The size of the system is given in n (n>=1).
  1188.  *
  1189.  * In the first procedure the Cholesky decomposition M = C^T * D * C
  1190.  * (C is upper triangle with unit diagonal, D is diagonal) is calculated.
  1191.  * Return TRUE if decomposition exist.
  1192.  */
  1193. static int solve_cubic_1(m, n)
  1194. tri_diag m[];
  1195. int n;
  1196. {
  1197.     int i;
  1198.     double m_ij, m_n, m_nn, d;
  1199.  
  1200.     if (n < 1)
  1201.     return FALSE;        /* Dimension should be at least 1 */
  1202.  
  1203.     d = m[0][1];        /* D_{0,0} = M_{0,0} */
  1204.     if (d <= 0.)
  1205.     return FALSE;        /* M (or D) should be positiv definite */
  1206.     m_n = m[0][0];        /*  M_{0,n-1}  */
  1207.     m_nn = m[n - 1][1];        /* M_{n-1,n-1} */
  1208.     for (i = 0; i < n - 2; i++) {
  1209.     m_ij = m[i][2];        /*  M_{i,1}  */
  1210.     m[i][2] = m_ij / d;    /* C_{i,i+1} */
  1211.     m[i][0] = m_n / d;    /* C_{i,n-1} */
  1212.     m_nn -= m[i][0] * m_n;    /* to get C_{n-1,n-1} */
  1213.     m_n = -m[i][2] * m_n;    /* to get C_{i+1,n-1} */
  1214.     d = m[i + 1][1] - m[i][2] * m_ij;    /* D_{i+1,i+1} */
  1215.     if (d <= 0.)
  1216.         return FALSE;    /* Elements of D should be positiv */
  1217.     m[i + 1][1] = d;
  1218.     }
  1219.     if (n >= 2) {        /* Complete last column */
  1220.     m_n += m[n - 2][2];    /* add M_{n-2,n-1} */
  1221.     m[n - 2][0] = m_n / d;    /* C_{n-2,n-1} */
  1222.     m[n - 1][1] = d = m_nn - m[n - 2][0] * m_n;    /* D_{n-1,n-1} */
  1223.     if (d <= 0.)
  1224.         return FALSE;
  1225.     }
  1226.     return TRUE;
  1227. }
  1228.  
  1229. /*
  1230.  * The second procedure solves the linear system, with the Choleky
  1231.  * decomposition calculated above (in m[][]) and the right side b given
  1232.  * in x[]. The solution x overwrites the right side in x[].
  1233.  */
  1234. static void solve_cubic_2(m, x, n)
  1235. tri_diag m[];
  1236. double x[];
  1237. int n;
  1238. {
  1239.     int i;
  1240.     double x_n;
  1241.  
  1242.     /* Division by transpose of C : b = C^{-T} * b */
  1243.     x_n = x[n - 1];
  1244.     for (i = 0; i < n - 2; i++) {
  1245.     x[i + 1] -= m[i][2] * x[i];    /* C_{i,i+1} * x_{i} */
  1246.     x_n -= m[i][0] * x[i];    /* C_{i,n-1} * x_{i} */
  1247.     }
  1248.     if (n >= 2)
  1249.     x[n - 1] = x_n - m[n - 2][0] * x[n - 2];    /* C_{n-2,n-1} * x_{n-1} */
  1250.  
  1251.     /* Division by D: b = D^{-1} * b */
  1252.     for (i = 0; i < n; i++)
  1253.     x[i] /= m[i][1];
  1254.  
  1255.     /* Division by C: b = C^{-1} * b */
  1256.     x_n = x[n - 1];
  1257.     if (n >= 2)
  1258.     x[n - 2] -= m[n - 2][0] * x_n;    /* C_{n-2,n-1} * x_{n-1} */
  1259.     for (i = n - 3; i >= 0; i--) {
  1260.     /*      C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1} */
  1261.     x[i] -= m[i][2] * x[i + 1] + m[i][0] * x_n;
  1262.     }
  1263.     return;
  1264. }
  1265.  
  1266. /*
  1267.  * Solve tri diagonal linear system equation. The tri diagonal matrix is
  1268.  * defined via matrix M, right side is r, and solution X i.e. M * X = R.
  1269.  * Size of system given in n. Return TRUE if solution exist.
  1270.  */
  1271. /* not used any more in "contour.c", but in "spline.c" (21. Dec. 1995) ! */
  1272.  
  1273. int solve_tri_diag(m, r, x, n)
  1274. tri_diag m[];
  1275. double r[], x[];
  1276. int n;
  1277. {
  1278.     int i;
  1279.     double t;
  1280.  
  1281.     for (i = 1; i < n; i++) {    /* Eliminate element m[i][i-1] (lower diagonal). */
  1282.     if (m[i - 1][1] == 0)
  1283.         return FALSE;
  1284.     t = m[i][0] / m[i - 1][1];    /* Find ratio between the two lines. */
  1285. /*      m[i][0] = m[i][0] - m[i-1][1] * t; */
  1286. /* m[i][0] is not used any more (and set to 0 in the above line) */
  1287.     m[i][1] = m[i][1] - m[i - 1][2] * t;
  1288.     r[i] = r[i] - r[i - 1] * t;
  1289.     }
  1290.     /* Now do back subtitution - update the solution vector X: */
  1291.     if (m[n - 1][1] == 0)
  1292.     return FALSE;
  1293.     x[n - 1] = r[n - 1] / m[n - 1][1];    /* Find last element. */
  1294.     for (i = n - 2; i >= 0; i--) {
  1295.     if (m[i][1] == 0)
  1296.         return FALSE;
  1297.     x[i] = (r[i] - x[i + 1] * m[i][2]) / m[i][1];
  1298.     }
  1299.     return TRUE;
  1300. }
  1301.  
  1302. /*
  1303.  * Generate a Bspline curve defined by all the points given in linked list p:
  1304.  * Algorithm: using deBoor algorithm
  1305.  * Note: if Curvekind is OPEN_CONTOUR than Open end knot vector is assumed,
  1306.  *       else (CLOSED_CONTOUR) Float end knot vector is assumed.
  1307.  * It is assumed that num_of_points is at least 2, and order of Bspline is less
  1308.  * than num_of_points!
  1309.  */
  1310. static void gen_bspline_approx(p_cntr, num_of_points, order, contr_kind)
  1311. struct cntr_struct *p_cntr;
  1312. int num_of_points, order, contr_kind;
  1313. {
  1314.     int knot_index = 0, pts_count = 1;
  1315.     double dt, t, next_t, t_min, t_max, x, y;
  1316.     struct cntr_struct *pc_temp = p_cntr, *pc_tail = NULL;
  1317.  
  1318.     /* If the contour is Closed one we must update few things:
  1319.      * 1. Make the list temporary circular, so we can close the contour.
  1320.      * 2. Update num_of_points - increase it by "order-1" so contour will be
  1321.      *    closed. This will evaluate order more sections to close it!
  1322.      */
  1323.     if (contr_kind == CLOSED_CONTOUR) {
  1324.     pc_tail = p_cntr;
  1325.     while (pc_tail->next)
  1326.         pc_tail = pc_tail->next;    /* Find last point. */
  1327.  
  1328.     /* test if first and last point are equal */
  1329.     if (fuzzy_equal(pc_tail, p_cntr)) {
  1330.         /* Close contour list - make it circular. */
  1331.         pc_tail->next = p_cntr->next;
  1332.         num_of_points += order - 1;
  1333.     } else {
  1334.         pc_tail->next = p_cntr;
  1335.         num_of_points += order;
  1336.     }
  1337.     }
  1338.     /* Find first (t_min) and last (t_max) t value to eval: */
  1339.     t = t_min = fetch_knot(contr_kind, num_of_points, order, order);
  1340.     t_max = fetch_knot(contr_kind, num_of_points, order, num_of_points);
  1341.     next_t = t_min + 1.0;
  1342.     knot_index = order;
  1343.     dt = 1.0 / num_approx_pts;    /* Number of points per one section. */
  1344.  
  1345.  
  1346.     while (t < t_max) {
  1347.     if (t > next_t) {
  1348.         pc_temp = pc_temp->next;    /* Next order ctrl. pt. to blend. */
  1349.         knot_index++;
  1350.         next_t += 1.0;
  1351.     }
  1352.     eval_bspline(t, pc_temp, num_of_points, order, knot_index,
  1353.              contr_kind, &x, &y);    /* Next pt. */
  1354.     add_cntr_point(x, y);
  1355.     pts_count++;
  1356.     /* As we might have some real number round off problems we do      */
  1357.     /* the last point outside the loop                                 */
  1358.     if (pts_count == num_approx_pts * (num_of_points - order) + 1)
  1359.         break;
  1360.     t += dt;
  1361.     }
  1362.  
  1363.     /* Now do the last point */
  1364.     eval_bspline(t_max - EPSILON, pc_temp, num_of_points, order, knot_index,
  1365.          contr_kind, &x, &y);
  1366.     add_cntr_point(x, y);    /* Complete the contour. */
  1367.  
  1368.     if (contr_kind == CLOSED_CONTOUR)    /* Update list - un-circular it. */
  1369.     pc_tail->next = NULL;
  1370. }
  1371.  
  1372. /*
  1373.  * The routine to evaluate the B-spline value at point t using knot vector
  1374.  * from function fetch_knot(), and the control points p_cntr.
  1375.  * Returns (x, y) of approximated B-spline. Note that p_cntr points on the
  1376.  * first control point to blend with. The B-spline is of order order.
  1377.  */
  1378. static void eval_bspline(t, p_cntr, num_of_points, order, j, contr_kind, x, y)
  1379. double t;
  1380. struct cntr_struct *p_cntr;
  1381. int num_of_points, order, j, contr_kind;
  1382. double *x, *y;
  1383. {
  1384.     int i, p;
  1385.     double ti, tikp, *dx, *dy;    /* Copy p_cntr into it to make it faster. */
  1386.  
  1387.     dx = (double *)
  1388.     gp_alloc((unsigned long) (sizeof(double) * (order + j)), "contour b_spline");
  1389.     dy = (double *)
  1390.     gp_alloc((unsigned long) (sizeof(double) * (order + j)), "contour b_spline");
  1391.  
  1392.     /* Set the dx/dy - [0] iteration step, control points (p==0 iterat.): */
  1393.     for (i = j - order; i <= j; i++) {
  1394.     dx[i] = p_cntr->X;
  1395.     dy[i] = p_cntr->Y;
  1396.     p_cntr = p_cntr->next;
  1397.     }
  1398.  
  1399.     for (p = 1; p <= order; p++) {    /* Iteration (b-spline level) counter. */
  1400.     for (i = j; i >= j - order + p; i--) {    /* Control points indexing. */
  1401.         ti = fetch_knot(contr_kind, num_of_points, order, i);
  1402.         tikp = fetch_knot(contr_kind, num_of_points, order, i + order + 1 - p);
  1403.         if (ti == tikp) {    /* Should not be a problems but how knows... */
  1404.         } else {
  1405.         dx[i] = dx[i] * (t - ti) / (tikp - ti) +    /* Calculate x. */
  1406.             dx[i - 1] * (tikp - t) / (tikp - ti);
  1407.         dy[i] = dy[i] * (t - ti) / (tikp - ti) +    /* Calculate y. */
  1408.             dy[i - 1] * (tikp - t) / (tikp - ti);
  1409.         }
  1410.     }
  1411.     }
  1412.     *x = dx[j];
  1413.     *y = dy[j];
  1414.     free((char *) dx);
  1415.     free((char *) dy);
  1416. }
  1417.  
  1418. /*
  1419.  * Routine to get the i knot from uniform knot vector. The knot vector
  1420.  * might be float (Knot(i) = i) or open (where the first and last "order"
  1421.  * knots are equal). contr_kind determines knot kind - OPEN_CONTOUR means
  1422.  * open knot vector, and CLOSED_CONTOUR selects float knot vector.
  1423.  * Note the knot vector is not exist and this routine simulates it existance
  1424.  * Also note the indexes for the knot vector starts from 0.
  1425.  */
  1426. static double fetch_knot(contr_kind, num_of_points, order, i)
  1427. int contr_kind, num_of_points, order, i;
  1428. {
  1429.     switch (contr_kind) {
  1430.     case OPEN_CONTOUR:
  1431.     if (i <= order)
  1432.         return 0.0;
  1433.     else if (i <= num_of_points)
  1434.         return (double) (i - order);
  1435.     else
  1436.         return (double) (num_of_points - order);
  1437.     case CLOSED_CONTOUR:
  1438.     return (double) i;
  1439.     default:            /* Should never happen */
  1440.     return 1.0;
  1441.     }
  1442. #ifdef sequent
  1443.     return 1.0;            /* ???? */
  1444. #endif
  1445. }
  1446.